Reaction Front in an A + B — ► C Reaction-Subdiffusion Process 
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We study the reaction front for the process A + B — » G in which the reagents move subdiffusively. 
Our theoretical description is based on a fractional reaction-subdiffusion equation in which both the 
motion and the reaction terms are affected by the subdiffusive character of the process. We design 
numerical simulations to check our theoretical results, describing the simulations in some detail 
because the rules necessarily differ in important respects from those used in diffusive processes. 
Comparisons between theory and simulations are on the whole favorable, with the most difficult 
quantities to capture being those that involve very small numbers of particles. In particular, we 
analyze the total number of product particles, the width of the depletion zone, the production profile 
of product and its width, as well as the reactant concentrations at the center of the reaction zone, 
all as a function of time. We also analyze the shape of the product profile as a function of time, in 
particular its unusual behavior at the center of the reaction zone. 

PACS numbers: 82.33.-z, 82.40.-g, 02.50.Ey, 89.75.Da 
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I. INTRODUCTION 

It is very well known that diffusion-limited binary re- 
actions in low dimensions may lead to the spontaneous 
appearance of spatial order and spatial structures, and 
to associated "anomalous" rate laws for the global den- 
sities of the reacting species. For example, the reac- 
tions A + A — > C (some selected references of many 
in the literature are QHHSSEHIlSGilllll) 
and A + B — > C (again, some of many references 

are|llllllllIi,|IllllElllillipi2ll2ll23) 

under "normal" circumstances are described by second- 
order rate laws, whereas the asymptotic rate law for the 
former reaction is of apparent order (1 + 2/d) for dimen- 
sions d < 2, and for the mixed reaction it is of apparent 
order (1 + 4/d) for d < 4. The slow-down implied by 
the higher order is a consequence of the rapid deviation 
of the spatial distribution of reactants from a random 
distribution. This is in turn a consequence of the fact 
that diffusion is not an effective mixing mechanism in 
low dimensions. 

To design an experiment in a constrained geometry 
in order to measure these anomalies is not at all sim- 
ple, especially for the mixed reaction [2~H I25L [2(j. It 
is simpler for the A + A problem because a number 
of appropriate non-chemical species can be identified 
that essentially undergo the simplest annihilation reac- 
tion or variants thereof. Examples include exciton an- 
nihilation experiments in one-dimensional por es and in 
effectively one-dimensional polymer wires excited 
molecule naphthalene fusion and quenching experiments 
in one-dimensional pores and kink-antikink simula- 
tions in one dimension 28]. Experimental observations 
of the A + B anomalies instead generally involve reac- 
tion fronts. Early on Galfi and Racz [2 ^ | and later oth- 
ers [MimmiHlHIHIilil recognized that 
the kinetic anomalies in the homogeneous systems would 
be reflected in the evolution of reaction fronts. On the 



basis of scaling arguments, later made more rigorous, a 
number of exponents were deduced to characterize this 
evolution. The first experiments confirming these results 
were carried out with species A and B diffusing in a gel 
contained in a thin capillary, with one species initially 
occupying one side, the other o ccup ying the other side, 
and a sharp front between them [24J . More recent exper- 



iments have been carried out in gel-free systems (25, 26 
The evolution of the reactant fronts and of the product of 
the reaction in these experiments both reflect the kinetic 
anomalies. 

In this paper we extend the front analysis to subdif- 
fusive reactions. Subdiffusive motion is characterized by 
a mean square displacement that varies sublinearly with 
time, 



(r 2 (t)) 



2K^ 



r(i + 7) 



(i) 



with < 7 < 1. For ordinary diffusion 7=1, and 
K\ = D is the ordinary diffusion coefficient. We argue 
for the importance of this generalization on a number 
of grounds. First, there exists a huge literature on sys- 
tems that deviate from diffusive behavior and are instead 
characterized by m otion all the way from subdiffusive to 
superdiffusive |3il l4Cj . Subdiffusive motion is particu- 
larly important in the context of complex systems such 
as glassy and disordered materials, in which pathways 
are constrained for geometric or for energetic reasons. It 
is also particularly germane to the way in which exper- 
iments in low dimensions have to be carried out. Such 
experiments must avoid any active or convective or ad- 
vective mixing so as to ensure that any mixing is only a 
consequence of diffusion. To accomplish this usually re- 
quires the use of gel substrates and/or highly constrained 
geometries (the first gel-free experiments were carried out 
recently |25|, |26(). Under these circumstances it is not 
clear whether the motion of the species is actually dif- 
fusive, or if it is in fact subdiffusive. Indeed, a recent 
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detailed discussion on ways to extract accurate parame- 
ters and exponents from such experiments concludes that 
at least the experiments presented in that work, carried 
out in a gel, reflect subdiffusive rather than diffusive mo- 
tion H3. 

We recently solved the A + A reaction-subdiffusion 
problem in one dimension |42|. To solve this problem, we 
generalized methods first applied to the reaction-diffusion 
A+A problem. For the A+A — > A problem the method of 
intervals allows an exact formulation in terms of intervals 
on the line that are empty of A particles 
The distribution of intervals evolves linearly, and there- 
fore one can find an exact solution. In the reaction- 
diffusion problem the description involves a diffusion 
equation, while the reaction-subdiffusion problem in- 
volves a subdiffusion equation j4^|; both can be solved 
exactly. For the A + A — > C problem one uses instead 
the odd/even parity method 0, whereby one keeps 
track of the parity of the number of particles in an inter- 
val. The associated distribution again satisfies a linear 
diffusion or subdiffusion [4^1 equation. 

Before these exact methods were developed, it was 
customary to model these and other binary reactions 
by writing down a reaction-diffusion equation for each 
species in the reaction. Such equations typically contain 
a diffusion term and a binary reaction term, the latter 
being some truncated form of a two-particle distribution 
function. For instance, in the A + B reaction the typical 
reaction term simply involves the product of the local 
concentrations of reactants, — ka(r, t)b(r, t) 0. In the 
A+A reaction one has to be slightly more careful because 
in writing, for example, — fca 2 (r,t) one must be careful 
not to include spurious self- reaction contributions 0, Q . 
Once the A + A exact models were developed that did 
not require one to explicitly write a reaction term, it 
was possible to analyze the accuracy of the approximate 
truncations Also, it was not necessary to consider 

the generalization of the reaction term to the subdiffu- 
sive case since the exact methods could be generalized 
directly. 

The situation is more complicated for the A + B prob- 
lem, because no such exact formulations or solutions have 
been developed in this case. There is a large literature on 
the reaction-diffusion problem with different truncation 
schemes to represent the reaction term, but the literature 
on the reaction-subdiffusion problem is far more recent 
and relatively unsettled. In particular, at the current 
stage of development of this problem it is necessary to 
think about how to (approximately) model the reaction 
term. 

In Sec. [n] we present a discussion of the model to be 
used for the description of the A+B reaction-subdiffusion 
problem. Having arrived at a particular set of fractional 
equations. We apply a scaling theory to these equations 
akin to that of Galfi and Racz [2^| , but now for a subdif- 
fusive front. To support the theoretical conclusions, it is 
necessary to perform numerical simulations, which is not 
a trivial matter for a problem involving subdiffusion. In 



Sec. lIIII we discuss our Monte Carlo simulation methods. 
Section II VI is a compendium and comparison of numer- 
ical and theoretical results. Some closing comments are 
presented in Sec. [V] 



II. THE MODEL 

We start with a system of A particles on one side and 
B particles on the other of a sharp linear front, defined 
to lie perpendicular to the x axis. The particles diffuse 
and react with a given probability upon encounter. A 
standard mean-field model for the evolution of the con- 
centrations a(x, t) and b(x, t) of A and B particles along 
x is given by the reaction-diffusion equations 



m a(x,t) = D-, 
d d 2 

—b(x,t) = D-^b(x,t) - ka(x,t)b{x,t) , 



, -a(x,t) - ka(x,t)b(x,t) 



(2) 



dt 



where D is the diffusion coefficient assumed to be equal 
for the two species. The initial conditions are that 
a(x, t) = const = ao for x < and a(x, t) = for 
x > 0. Similarly, b(x, t) — const — bo for x > 
and b(x,t) — for x < 0. With these conditions, 
no matter the dimensionality of the system, the sys- 
tem of equations is effectively one-dimensional. The 
front problem was first analyzed via a scaling descrip- 
tion |29| and later refined by a large number of authors 
using more ri goro us theoretical and careful numerical ap- 
proaches [3(J, IU H2, yj, |35|, |3a, |33, |33| . One upshot 
of the extensive work is that d = 2 is a critical dimension 
for the mean field description to be appropriate. Below 
d — 2 one must take into account fluctuations, neglected 
in this description, that completely change the outcome 
of the analysis. A particularly transparent argument for 
this critical dimension was provided by Krapivsky [3^| . 
He argued that the reaction constant in the mean-field 
reaction rate r = kab should in general depend on the 
diffusion constant D and the radius R of the reacting 
particles. Dimensional analysis gives k ~ DR d ~ 2 , but 
on physical grounds one expects the reaction rate con- 
stant to be an increasing function of the radius R. The 
conclusion is that the mean field model can therefore not 
be valid for d < 2. While it has been assumed that the 
mean field model holds for the critical dimensions d = 2, 
Krapivsky finds logarithmic corrections that have also 
been observed in simulations [3]]]. In our analysis and 
simulations we will take d = 2 (which turns out to be the 
critical dimension for the subdiffusive problem as well) 
and will therefore not deal with the lower-dimensional 
fluctuation effects. In this first study we will not deal 
with logarithmic corrections. 

In order to generalize the reaction-diffusion problem to 
reaction-subdiffusion, we must deal with the subdiffusive 
motion of the particles (generalization of the first term in 
Eq. (J2J) and with their reaction rate law (second term). 
We discuss each separately. 
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SubdifFusion is not modeled in a universal way in the 
literature. Among the more successful approaches to 
the subdiffusion problem have been continuous time ran- 
dom walks with non-Poissonian waiting time distribu- 
tions 0, EH Ell , and fractional dynamics approaches in 
which the diffusion operator is replaced b y a g eneralized 
fractional diffusion operator [39LI45L l47ll4a l49 l \ The rela- 
tion between the two has also been discussed 39, 4511431. 



In particular, the fractional dynamics formulation that 
leads to the mean square displacement Q can be associ- 
ated with a continuous time random walk with a waiting 
time distribution between steps which at long times be- 
haves as 



(3) 



We adopt the fractional dynamics approach, and com- 
ment later on some issues associated with it that must 
carefully be considered in the context of numerical simu- 
lations. We thus replace Eq. J5J) with the set of reaction- 
subdiffusion equations 

d d 2 

—a(x,t) = K 1 dI~~< -^a(x,t) - R 1 {x,t) 

d d 2 

— b(x,t) = K 1 oD 1 ^ '—b(x,t) -Ry(x,t) , 

where K~ is the generalized diffusion coefficient that ap- 
pears in Eq. |JT||. and o^* -7 is the Riemann-Liouville 
operator, 



a 1_7 /(m) 



1 8 f\ T JM 



r(7) dt Jo (t - tY~i 



(5) 



The reaction term _R 7 (x, t) = R will be discussed subse- 
quently, because certain aspects of the problem are inde- 
pendent of the specific form of this term. 



A. Scaling independent of reaction term 

As the reaction proceeds, a depletion zone develops 
around the front. This is the region where the concen- 
trations of reactants are significantly smaller than their 
initial values. How the width Wd evolves with time is 
one of the measures typically used to characterize the 
process. Within this depletion zone lies the so-called re- 
action zone, the region where the concentration c(x, t) of 
the product C is appreciable. This concentration profile 
has a width w whose variation with time is another char- 
acteristic of the evolving reaction. The evolution of the 
production rate of C (which determines the height of the 
profile of c(x, t) in the reaction zone) is a third measure 
of the process. To find these time dep endences we adapt 
the original scaling approach [29l l3l| to the subdiffusivc 
case, and assume the scaling forms 



a(x,t) = t- 9 a(xt~ a ) 
b(x,t) = t' b(xt- a ) 



(6) 



for the concentrations and 

Ry(x,t) = t-»Ry(xr 



(7) 



for the reaction term. The exponents 9, a, and /x are to 
be determined from three relations. The scaling forms are 
only valid for x <C Wd, that is, well within the depletion 
zone. 

Two of the three relations needed to fix the scaling ex- 
ponents do not require further specification of the reac- 
tion term. Since the reaction zone increases more slowly 
than the width of the depletion zone (an assumption that 
ex post turns out to be correct), we can focus on the 
concentration difference u(x,t) = a(x,t) — b(x,t) to de- 
duce the width of the latter. The reaction term drops out 
when one subtracts the equations in Eq. J2J , and its form 
therefore does not matter at this point. Generalizing the 
procedure in [2^| , one can scale the resulting equation by 
measuring concentrations in units of ao, time in units of 



t = l/(fcao), and length in units of / 
that the equation is simply 



(if 7 T 7 ) 1/2 : SO 



d_ 

Of 



u(x, t) 



dx 2 



(8) 



and the only control parameter is q = bo/ao in the initial 
condition: 

u(x, 0) = 1 for x < 
u(x, 0) = — q for x > 0. 



(9) 



The solution is 



u(x, t) 



q n i,o 



X 


(i.i)" 


rii 2 


(0,1) _ 



where H 1 \ is the Fox H-function |4' 
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7=1 this reduces to the diffusion result [2^| 

i + g 



(10) 



When 



u{x, i) 



-Q 



erfc 



(2^/2) 



(11) 



From Eq. i|10|l we see that the width of the depletion zone 
scales as 



W d ~ v' 2 



(12) 



i.e., da(x,t)/dx ~ db(x,t)/dx ~ V 1 ! 2 . Then, from 
Eq. © , the following relation between scaling exponents 
follows immediately: 



7 

a = — . 
2 



(13) 



The second relation follows from the fact that the con- 
centration gradient of A and B leads to a flux of particles 
toward the reaction region. The assumption that the re- 
action is fed by these particle currents then leads to the 
quasistationary form in the reaction zone 



= K 7 oDj-fpo^t) 
= K 7 D]-^b(x,t) 



- R 1 {x, t) 



(14) 
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which requires that 

^ = 6» + 2a + l- 7 



(15) 



For the width of the reaction zone to grow more slowly 
than the depletion zone caused by subdiffusion requires 
that 



<7/2 



(16) 



On the other hand, the quasistationarity condition re- 
quires that 



dx 2 



» jLa(x,t) ~ t-^ +1 \ (17) 



which again leads to Eq. i|16f) . 

Equations l|13|) and l|15|) combined lead to the relation 
a — [i = 7/2 — I that is easily checked by numerical 
simulations, since it is determined by the production rate 
of C . The rate of change of the total amount of product, 
dNc/dt, is given by the integral of the reaction rate over 
the reaction zone, 

dNn 



, dx Ry(x,t) ~ / dxRj{x/t a ) 

dt J reaction / reaction 



t-0-«) ~fr/2-i ; (18) 



that is, 



N c (t) ~ i 7/2 



(19) 



We stress that this total amount of product as a func- 
tion of time, which is numerically more robust than its 
derivative, is predicted to grow as t 1 ! 2 regardless of the 
specific form of the reaction term. 

Another accessible quantity that is independent of i? 7 
is the location Xf of the point at which the production 
rate of C is largest. This should occur where a(x,t) ~ 
b(x, t), that is, u(xf, t) ~ 0. The time dependence of this 
equimolar point is found from Eq. @ to be 



x f (t)=K f t^ 2 
where Kf is determined from the equation 



_2q_ 
l + q 



ttW 

H ll 



Kf 



(Li) 
(0,1) 



(20) 



(21) 



memory with the reaction term. Some assume that, as 
in the case of ordinary diffusion, reactions can simply be 
modeled by a space-dependent form of the law of mass 
action, e.g., by setting R = ka(x, t)b(x, t). Some of these 
assumptions may be appropriate if the reaction is very 
rapid, but not if many encounters between reactants are 
required for the reaction to occur. 

We adopt the viewpoint put forth in a recent theory 
developed for geminate recombination 0, |5l| but, as 
the authors themselves point out, much more broadly 
applicable. This theory goes back to the continuous time 
random walk picture from which the fractional diffusion 
equation can be obtained, and considers both the mo- 
tion and the reaction in this framework. In the context 
of geminate recombination the authors define a reaction 
zone and argue that a geminate pair within the reaction 
zone will not necessarily react for any finite intrinsic re- 
action rate (which they call j rc ) because one of the par- 
ticles may leave the zone before a reaction takes place. 
The dynamics of leaving the reaction zone is ruled by 
the waiting time distribution ^p ut{t) = i/'(0 e_7rc ' where 
ip(t) is the waiting time that regulates the rest of the dy- 
namics [cf. Eq. , and therefore the reaction rate will 
acquire a memory that arises from the same source as 
the memory associated with the sub-diffusive motion. In 
the continuum limit this model then leads to a reaction- 
subdiffusion equation in which both contributions have 
a memory. Seki et al. obtain a subdiffusion-reaction 
equation which at long times corresponds to choosing a 
reaction term of the form 



R 1 {x,t)=kaD t 1 a(x,t)b(x,t). 



(22) 



Here "long times" set in very quickly if the reaction zone 
is narrow and the intrinsic reaction rate small. As noted 
earlier, although the derivation is specifically for gemi- 
nate recombination, the arguments can be generalized. 

Our full reaction-subdiffusion starting equations on 
which the remainder of this paper is based then are 



d ( d 2 

—a(x,t) = oA -7 \ K ig^ a ( x ' t ) ~ ka(x,t)b(x,t) 



d 



^-b(x,t) = dI 7 i K 7 -^b(x,t) ka{x,t)b(x,t) 



dt 



o 2 



dx 2 



(23) 

From the specific reaction term given in Eq. (12211 we can 
now obtain the third relation between the scaling expo- 
nents by balancing the terms within the brackets: 



B. Choice of reaction term and resultant scaling 



H = 26 + 1 - 7. 



(24) 



Further relations involving the scaling exponents 
aimed at their expression in terms of model quantities re- 
quire specification of the reaction term. There is a varied 
literature on this subject, based on a number of different 
assumptions Ejjl 03 E3- Most do not associate a 



Simultaneous solution of Eqs. I|13fl . (|15fl . and l|24|l finally 
yields 
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C. Simulated quantities 



III. SIMULATION DETAILS 



It is useful to list here the quantities that will be com- 
pared with numerical simulations. Each is character- 
ized by an exponent explicitly given in terms of 7. The 
first and second are independent of the choice of reaction 
term, but the others are sensitive to this choice. 

1 . The total amount of product C produced as a func- 
tion of time, given in Eq. (|19|l . is 



N c (t) ~ t^ 2 



(26) 



This scaling is independent of the form of the reac- 
tion term. 

2. We measure the width Wd of the depletion zone as 
the width of the profile 



U P (x,t) = 1 - \a(x,t) - b(x,t)\. 



(27) 



The prediction, which is also independent of the 
form of the reaction term, is given in Eq. (|12|l . 



W d ~ V' 2 



(28) 



3. We carry out our simulations with an equal initial 
unit concentration of A and B. In this case Xf = 
for all time. We monitor the number of C parti- 
cles produced at this point of maximum produc- 
tion of C, N c {x = 0, t). Since i? 7 (0, t) = dN c {x = 
0,t)/dt ~ t-f* = p/ 3 " 1 , we have 



N c (0,t) ~F /3 . 
This is thus a check on the exponent \x. 



(29) 



4. The concentration a(0, t) = 6(0, t) of each reactant 
at the center of the reaction zone is difficult to mon- 
itor because it is very small and therefore subject to 
large fluctuations. Instead, we monitor the integral 
of this concentration over time, 



t rt 
11 Jo 



(30) 



with a similar result for the other reactant. This 
then is a check on the exponent 9. 

5. The width w(t) of the product profile grows, ac- 
cording to the scaling equation Q, as w(t) ~ t a . 
According to Eq. (|25|l we then have, as a test of a, 



■;(t) ~ T /6 



(31) 



6. Finally, we monitor the entire profile (|27|l as a func- 
tion of position and time. This is a difficult quan- 
tity to follow because it involves regions of very low 
concentration. In a way it constitutes a check of the 
simulation methodology, as we will see below. 



Monte Carlo simulation methods for reaction-diffusion 
processes are ubiquitous. For a two-dimensional simula- 
tion one starts with a square lattice, and deploys a given 
number of particles at each site according to the initial 
distribution. The particles then perform a random walk 
simulated by the parallel update of the coordinates of all 
particles at each time step t — mAt, m = 1, 2, • • • . The 
entire lattice is explored at periodic intervals At r (which 
could and often does coincide with At), and reactions 
take place at each site on which there are A and B par- 
ticles, with probability kab. Here k <C 1 is the reaction 
rate constant and a and b are proportional to the number 
of particles of type A and B on that site. Clearly, kab 
must (in the appropriate sense, since the quantity is not 
dimensionless) be small. There are variants of this proce- 
dure that are inconsequential for our analysis (e.g., some 
excluded volume effects). A necessary condition to be 
in the diffusion-controlled regime described by the usual 
reaction-diffusion equations is that the random walkers 
on average perform a large number of steps before react- 
ing. 

Adjustments that must be made to this procedure 
in order to describe subdiffusion are neither trivial nor 
straightforward. First and most importantly, one can not 
assume that the particles all jump at the same time. The 
distribution of jumping times is now very broad: one can 
imagine each particle outfitted with an alarm clock, with 
a jump to a randomly selected nearest neighbor taking 
place when the alarm goes off, at which time the alarm 
is reset according to a distribution whose asymptotic be- 
havior goes as in Eq. Jumping is therefore a renewal 
process [Hfij . An example of a normalized distribution 
with this behavior is the Pareto law: 



7/r 



(32) 



The particles are labeled, and jumping times are assigned 
to them according to this distribution. These times, from 
smallest to largest, must be sorted, and the list must be 
sorted after each jump or reaction. 

Since the particles no longer jump at the same time, a 
decision must be made about when they are allowed to 
react. There are at least two alternatives: (1) A reac- 
tion attempt occurs only when a particle first arrives at 
a site; (2) Reaction attempts occur at each site at peri- 
odic intervals At r , and occur with probability kab. The 
first alternative does not seem physically reasonable for 
the subdiffusive problem since it implies that a pair of 
A and B particles that remain at a given site and that 
did not react upon first encounter will not react no mat- 
ter how long they remain at the site, which on average 
is infinite. They can only react if they move apart and 
then encounter one another again. The second alterna- 
tive, which we choose for most of our simulations, can be 
associated with a number of physical explanations. On 
the one hand one can think of reactions induced or ac- 
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tivated periodically by some external agent (a laser, for 
example). More in line with our thinking of subdiffusion 
as a way to describe movement in a disordered or glassy 
or porous medium is to think of this as a mesoscopic de- 
scription. At a microscopic level small jumps may occur 
diffusively, but the motion from one mesoscopic region to 
another on a longer time scale is much slower because of 
geometric bottlenecks that affect this longer range move- 
ment. Our "sites" would then correspond to mesoscopic 
regions in which a walker can spend a long time moving 
diffusively from one part of the region to another. Reac- 
tions can then take place within one of these regions at 
regular time intervals. 

Since the subdiffusive process has a long memory, we 
must be careful about the initiation of the process. In 
particular, it is not appropriate to choose the initial 
jumping times as indicated above because that would 
bias the initial condition to one in which all the particles 
jumped simultaneously at time t = 0. Instead, after this 
first selection of times we choose another set of jumping 
times from the distribution, and repeat this procedure 
a large number of times. The number of repetitions is 
usually chosen as the total number of particles initially 
in the system. Only then do we choose t = by taking 
the smallest jumping time as our new origin of time from 
which the process is launched. 

Finally, it is noted that even in a diffusive process, 
reaction events are not really restricted to occur only at 
periodic time intervals At r . A large literature points to 
the fact that the continuous time process underlying such 
a step process is one in which times are selected from an 
exponential distribution |57ll58| 

P(t) = Ke- K \ (33) 

where re is the reaction rate constant. We have also tested 
this procedure in our subdiffusive system, allowing each 
pair (A, B) of particles on one site to react at a time dic- 
tated by such an exponential distribution. If a particle 
leaves a site before a reaction takes place, the reaction 
"clock" of each particle is reset. This is also the view- 
point followed by Seki et al. [3, EH- Note that whereas 
in the k, At r formulation of the reaction events one spec- 
ifies two parameters, in the exponentially distributed re- 
action events the reactions are specified by the single rate 
parameter re, which identifies both the reaction rate con- 
stant re and the mean time between reaction events 1/re. 

The parameters used in our simulations are as follows: 
we place two particles initially at each site (this corre- 
sponds to unit concentration for each species). The lat- 
tice dimensions are usually (L x ,L y ) = (20,10), except 
for 7 = 3/4 where we use L x = 90, and in some cases 
specified later where we use (L x ,L y ) = (160,20). The 
maximum number of particles allowed at a given site is 
40. The rate coefficient is k = 0.05, and the time be- 
tween reaction events is At r = 10. For some of the sim- 
ulations we use exponentially distributed reaction events 
with re = 10~ 4 or re = 10~ 5 , which corresponds to a 
much lower reaction rate. The maximum time per run is 
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FIG. 1: Log-log plots of the total number of product particles 
vs time for 7 = 0.75 (squares), 7 = 0.5 (circles), and 7 = 
0.25 (triangles). The linear fit slopes are 0.37, 0.25, and 0.15 
respectively. The mean field prediction for the slope is given 
in Eq. as 7/2. 
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FIG. 2: Width Wd of the depletion zone vs time for 7 = 0.5 
(circles) and 7 = 0.25 (triangles). The linear fit slopes are 
0.257 and 0.154 respectively. The mean field prediction for 
the slope is given in Eq. 1281 as 7/2. 

tmax — 1,024,000. Results are averaged over 100 runs. 

IV. COMPARISONS WITH SIMULATIONS 

Here we compare simulations of the six quantities 
enumerated in Sec. [H] with the theoretical predictions. 
The simulations rapidly become increasingly difficult and 
time-intensive with decreasing 7, and it is therefore ex- 
pected that agreement with the theory improves with 
increasing 7. As we will show, the agreement is on the 
whole good, especially for the larger values of 7. We also 
stress that four of the six comparisons involve results that 
decidedly depend on the choice of reaction term. Agree- 
ment would not be obtained with the usual memoryless 
local law of mass action. 

Figure^shows our simulation results for the total num- 
ber of product particles as a function of time in units of 
r = 1.0 (used throughout) for several values of 7. This is 
perhaps the most robust global quantity to be simulated. 
The linear fit slope is given for each 7, and agrees very 
well with the theoretical prediction given in Eq. (|26|l for 
the two larger values of 7. 

Figure [21 shows our simulation results for the width of 
the depletion zone as a function for time for two values 
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FIG. 3: Log-log plots of the production profile of product C 
at x = as a function of time for 7 = 0.75 (squares), 7 = 0.5 
(circles), and 7 = 0.25 (triangles). The linear fit slopes are 
0.24, 0.162, and 0.093 respectively. The mean field prediction 
for the slope is given in Eq. 129H as 7/3. 
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FIG. 4: Time integrals of reactant concentrations at the cen- 
ter of the reaction zone. The mean field prediction for the 
slopes is given in Eq. I)30|l as 1 — 7/3. The steeper curve is for 
7 = 0.25 and the reaction events governed by the exponen- 
tial distribution Eq (1331 with k = 10~ 5 . The linear fit slope 
is 0.814 while mean field theory yields 0.917. The shallower 
curve is for 7 = 0.5 and k — 10 -4 . The linear fit slope is 
0.683, the mean field slope is 0.833. 
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FIG. 5: Log-log plot of the width of the product profile as a 
function of time for 7 = 0.75 (squares), 7 = 0.5 (circles), and 
7 = 0.25 (triangles). The linear fit slopes are 0.129, 0.084, and 
0.042 respectively. The mean field prediction for the slope is 
given m Eq. as 7/6. 
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FIG. 6: Simulation results for the profile Up(x, t) for 7 = 0.75 
and t = 946, 3777, 15073, 60149, 240025, and 957828. The 
width of the profiles increase with time, as seen in Fig. 
The simulation was carried out on a lattice of size (160,20) 
and averaged over 102 runs. Notice the apparent evolution 
from a characteristic sharp-pointed profile for short times to 
a vaulted profile at longer times. For a discussion of this 
anomaly, and for the values of other parameters, see text. 



sion of the profile Up(x,t). 



of 7. The linear fit slope is in good agreement with the 
theory as given in Eq. I|28|) for the larger value of 7. Later 
we discuss some difficulties, particularly for small values 
of 7, in the accurate simulation of the profile Up(x,t) 
whose width is used to obtain these results. 

Figure |3 shows our simulation results for the produc- 
tion profile of the product of the reaction at x — as a 
function of time, for several values of 7. For the larger 7 
the linear fit slope agrees well with the theoretical pre- 
diction given in Eq. I|29|l . 

Figure 0] presents our simulation results for the time- 
integral of the concentration of a reactant at the center 
of the reaction zone, to be compared with the theoret- 
ical prediction of Eq. H3U|) . While the agreement is not 
spectacular, the trend is correct. Also of interest here is 
the improved agreement when the reaction rate is greatly 
reduced, as expected. Even more dramatic effects of the 
reaction rate are seen below in the context of our discus- 



Figure |S] contains our simulation results for the width 
of the product profile, which should be compared with 
the prediction of Eq. (|31l) . The agreement is very good 
for all the values of 7. 

Finally, in Figs. El and [7] we present perhaps the most 
difficult quantity to capture accurately, namely, the pro- 
file U p (x, t) defined in Eq. I|27|) . For these simulations we 
use a lattice of size (L x , L y ) = (160, 20). As pointed out 
earlier, the difficulty arises from the fact that it involves 
regions of very low concentration. It is instructive to il- 
lustrate the difficulty, and that is why we have included 
Fig. HJ1 The simulation profiles shown at the different 
times are actually time averaged over a small time inter- 
val around the times shown. The noteworthy feature is 
the evolution of the sharp-pointed profile near the origin 
at short times to a more rounded shape at longer times. 
The mean field theory presented in this paper does not 
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FIG. 7: Simulation (jagged lines) and theory (dashed lines) 
for the profile U P (x, t) for 7 = 0.75 and t = 3738, 8577, 19682, 
45165, 103638, 237815, and 545708. Notice the absence of the 
rounding anomalies in the profile for small x. The reaction 
rate is much smaller in this case, see text. The simulation was 
carried out on a lattice of size (160, 20) and averaged over 102 
runs. The times are approximately equivalent to those used 
in Fig. [fj] if the total number of C particles in the system is 
used as a measure of time. 

produce this rounding, so that it seemed at first that the 
theory and simulations differed in some profound way. 
However, the simulations in Fig. H3 were carried out with 
the rate coefficient k = 0.05 with reactions occurring pe- 
riodically at time interval At r = 10, a reaction rate that 
turns out to be too high for comparison with our theory. 
In Fig. [7] we show both the simulation results (jagged 
curves) and those of our theory (dashed curves), now 
with exponentially distributed reaction events according 
to Eq. with k = 10 -5 . The profile near x = 

remains pointed for all the times shown, as predicted 
by the theory. The quantitative disagreements for long 
times (t = 237185 and especially t = 545708) are due to 
finite size effects. Boundary effects are negligible only as 
long as U p (x = —L x /2,t) = Up(x = L x /2,t) = 0, and 
for the long times we find that this is not the case. In 
our other simulation results we have not included such 
results in our averages, but have left them in this fig- 
ure simply to stress the finite size effects one must be 
aware of in calculating the behavior of quantities in the 
depletion zone when the zone extends all the way to the 
boundaries of the system. To provide accurate results 
for such long times it is necessary to run simulations on 
larger systems. 

V. CLOSING COMMENTS 

In this paper we have proposed a set of continuum 
fractional diffusion equations to describe the behavior of 
a reaction front in the A + B — » C reaction- subdiffusion 
problem. Subdiffusion may be appropriate to describe 
the way reactants move in complex (glassy, disordered, 
highly constrained) geometries, and we were interested in 
exploring how this constraint on the motion would affect 



the evolution of a reaction front. Because we are working 
with a set of mean field continuum equations, our results 
are only valid above the critical dimension d = 2. 

The subdiffusive motion is modeled via the usual frac- 
tional equation that contains the Riemann-Liouville op- 
erator, Eq. J5J. This choice has a long history, and its 
virtues and shortcomings are clearly understood. Less 
clear has been the selection of the local reaction term, 
and the question of the way in which the memory in 
the Riemann-Liouville operator does (or does not) af- 
fect the way in which the reaction is modeled. While 
the literature on this subject has presented a number of 
viewpoints, we argued, in agreement with |54l l55| . that 
the reaction term should also be modified from its usual 
simple instantaneous product form, at least for small re- 
action rate constants. Our reaction-subdiffusion model 
is thus given by Eq. l(2"3l . 

Following the approach of Galfi and Racz [2^| for the 
evolution of a front in the reaction-diffusion problem, we 
assumed scaling solutions for the various quantities that 
can be calculated from the model. Some of these quanti- 
ties depend explicitly on the form chosen for the reaction 
term while others do not. We compared the resulting ex- 
ponents with those obtained from numerical simulations. 
We found very good agreement between the theory and 
simulations for the exponents /1 and a that characterize 
the reaction term, Eq. J7J. In particular, in terms of the 
power 7 < 1 that characterizes the subdiffusive process 
we found that the theoretical values [i = 1 — 7/3 and 
a = 7/6 are recovered in the simulations with greater 
fidelity for larger 7. The exponent 9, governing the time 
decay of the reactant concentrations as in Eq. JBJ, is, 
theoretically, given by 7/3. Simulation results give val- 
ues which correctly follow this trend, but the agreement is 
not quantitative. However, we have to remember that the 
quantity a(0, t) is local and, consequently, it is more dif- 
ficult to achieve good statistical averages with the small 
systems and number of particles considered in the sim- 
ulations. Perhaps the most challenging quantity to cap- 
ture is the profile Up(x,t). The theory predicts a cusp 
at x = which we were able to capture by our simula- 
tions when the reaction rate constant is sufficiently small. 
The quantitative agreement between the theory and the 
simulations for this profile was ultimately limited by our 
finite system size. We note that our results are a good ex- 
ample of what is sometimes referred to as subordination 
in that the subdiffusive scaling behavior can be deduced 
from the corresponding diffusive behavior with the sub- 
stitution t-fiT |59|. 

This work can clearly be pursued along a number of 
directions. Among them is the description of this same 
reaction-subdiffusion problem with the usual uniform ini- 
tial condition for the species, to investigate what sorts 
of segregation patterns might evolve on the way to ex- 
tinction, or on the way to equilibration if the reaction 
is reversible. Another is the study of the fluctuations 
that must be added to the model in order to describe 
this process in a one-dimensional system where the mean 
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field description is no longer appropriate, and the pos- 
sible logarithmic correction in two dimensions that may 
explain some of our small-7 discrepancies. A third is the 
effect of different subdiffusion coefficients for the species 
A and B, and even of different exponents ja and jb- In 
this latter case subordination would necessarily be more 
complicated if valid at all. Work along these directions 
is in progress 60] . 
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